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Abstract 

This paper describes a multidimensional hydrodynamic code which can be used 
for studies of relativistic astrophysical flows. The code solves the special relativistic 
hydrodynamic equations as a hyperbolic system of conservation laws based on the 
total variation diminishing (TVD) scheme. It uses a new set of conserved quantities 
and employs an analytic formula for transformation from the conserved quantities 
in the reference frame to the physical quantities in the local rest frame. Several 
standard tests, including relativistic shock tubes, a relativistic wall shock, and a 
relativistic blast wave, are presented to demonstrate that the code captures dis- 
continuities correctly and sharply in ultrarelativistic regimes. The robustness and 
flexibility of the code are demonstrated through test simulations of the relativistic 
Hawley-Zabusky shock and a relativistic extragalactic jet. 
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1 Introduction 



Many high-energy astrophysical problems involve relativistic flows, and thus 
understanding relativistic flows is important for correctly interpreting astro- 
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physical phenomena. For instance, intrinsic beam velocities larger than 0.9c 
are typically required to explain the apparent su perluminal motions observed 
in relativistic jets in microquasars in the Galaxy ([Mirabel fc Rodrigue 
as well as in ex tragalactic radio sources associated with active galactic nuclei 
( ZensusL Il997h . In some powerful extragalactic radio sources, ejections from 
galactic nuclei produce true beam velocities of more than 0.98c. Relativistic 
descriptions are also inevitable in other s ituations of rap id expansion such as 
the early stages of supernova e xplosions (iBurrowsj . ioOO) and the production 
of energetic gamma-ray bursts f M eszarosl . 120021 ) . In the general fireball model 
of gamma-ray bursts, the internal energy of gas is converted into the bulk ki- 
netic energy during expansion and this expansion leads to relativistic outflows 
with high bulk Lorentz factors > 100. Since such relativistic flows are highly 
nonlinear and intrinsically complex, in addition to possessing large Lorentz 
factors, often studying them numerically is the only possible approach. 



For numerical study of non-relativistic hydrodynamics, explicit finite differ- 
ence upwind schemes have been developed and implemented successfully. The 
schemes which hay e been used for astrophysical researches include the Roe 
schem e ( RoeL 1981 ). the total variation diminishing ( TVD) scheme (lliarten . 

1983 ) . the piecewise parabolic method (PPM) scheme (IColella k. Woodward . 

1984 ) . and the essentially non-oscillatory (ENO) scheme fllarten et all Il987l ) 
These schemes are based on exact or approximate Riemann solvers using the 
characteristic decomposition of the hyperbolic system of hydrodynamic con- 
servation equations. They all are able to capture sharp discontinuities robustly 
in the complex flows, and to describe the physical solution accurately. 



Although the upwind schemes were originally developed for non-relativistic 
hydrodynamics, s ome have been ex t ended to special relativistic hydrodynam- 
ics. For instance, iDolezal fc Wong (1995) adapted the ENO scheme to one- 
dimensional relativistic hydrodynamics. They fulfilled the ENO scheme using 
the local characteristic approach w hich depends on t he loc al linearizion of the 
system of conservation equations. Marti fc Miiller (1996) adapted the PPM 
scheme to one-dimensional relativistic hydrodynamics using an ex act relativis- 



tic Riemann solver to calculate numerical fluxes at cell interfaces. Donat et al. 



( 19981) and lAlov et all (|l999l) constructed multidimensional relativistic hydro 
dynamic codes based on the ENO scheme and the PPM scheme, respectively. 
Reviews of various nu meric al approaches and test pro blems can be found in 
Marti fc Mulierl (|2003h and IWilson fc Mathewsl (|2003h . These works showed 
that the advantage of the upwind schemes, high accuracy and robustness, are 
carried over to relativistic hydrodynamics. 



In this paper we describe a multidimensional code for special relativ istic hy- 
drody namics based on the total variation diminishing (TVD) scheme (IHarten . 
1983). The TVD scheme is an explicit Eulerian finite difference upwind scheme 
and an extension of the Roe scheme to second-order accuracy in space and 
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time. The advantage of the TVD scheme is that a code based on it is simple 
and fast, and yet performs well. A non-relativistic hydrodynamic code based 
the TVD scheme was built and applied to astrophysical problems such as the 



large scale structure formation of the universe by one of authors (|Rvu et al 



1993). The special relativistic hydrodynamic code in this paper was built 



by extending the non-relativistic code. All the components of the the non- 
relativistic code was kept, so the relativistic code has the structure parallel to 
the non-relativistic counterpart. It makes the relativistic code comprehensible 
and easily usable. Through tests, we demonstrate that the newly developed 
code for special relativistic hydrodynamics can handle interesting astrophysi- 
cal problems involving large Lorentz factors or ultrarelativistic regimes where 
energy densities greatly exceed rest mass densities. 

This paper is organized as follows. In Section 2 we describe the step by step 
procedures for building the code including the basic equations, characteristic 
decomposition, TVD scheme, multidimensional extension, and Lorentz trans- 
formation. Tests are presented in Section 3. A summary and discussion follows 
in Section 4. 



2 Numerical Relativistic Hydrodynamics 



2.1 Basic Equations 



The ideal relativistic hydrodynamic equations can be written as a hyperbolic 
system of conservation equations 

dM- d 

^ + ^-(M^,+^.) = 0, (2) 
_ + [(*+,)„] = 0, (3, 



where the equation of state is given by 

p=( 7 -l)(e-p). (4) 



Here, D, M i} and E are the mass density, momentum density, and total energy 
density in the reference frame, and p, Vj, and e are the mass density, velocity, 
and internal plus mass energy density in the local rest frame, respectively. In 
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general, the adiabatic index 7 is taken as 5/3 for mildly relativistic cases and 
as 4/3 for ultrarelativistic cases where e ^> p. In equations (l)-(3), the indices 
i and j run over x, y, and z and the conventional Einstein summation is used. 
The speed of light is set to unity (c = 1) throughout this paper. 

The quantities in the reference frame are related to those in the local rest 
frame via Lorentz transformation 



D = Tp, 

Mi = T 2 (e + p) Vi , 
E = T 2 (e + p)-p, 



(5) 
(6) 
(7) 



where the Lorentz factor is given by 
1 



(8) 



with v 2 = v\ + v 2 + v\. 



In the non-relativistic limit, the quantities D, Mi, and E approach their non- 
relativistic counterparts p N , p N vf, and E N + p N c 2 and equations (l)-(3) re- 
duce to the non-relativistic hydrodynamic equations 



dp 



N 



d 



dt dx 



d P N v? 

dt 
dE N 

+ 



d 



dt 



dx 
d 

dxj 



E N +p N )vf 



0, 



(9) 
(10) 

(11) 



where the pressure is given by 

P N - (7 - 1) (> - \p N v N2 ) • 



'12) 



2.2 Characteristic Decomposition 



Equations (l)-(3) can be written as 



dq dFj 
— H = 

dt dxi 



(13) 
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with the state and flux vectors 





' D ' 




Dv 3 


Q = 


Mi 


, F S = 


MiVj + pSij 




E 




. (E + p) Vj _ 



(14) 



or as 



dt ^ 3 dxi 



A, = -A 



(15) 



Here, A, is the 5x5 Jacobian matrix composed with the state and flux vec- 
tors. The construction of the matrix Aj can be simplified by introducing a 
parameter vector, u, as 



A A 



dFj du 
du dq 



(16) 



We choose the parameter vector which consists of the physical quantities in 
the local rest frame, 



u 



P 

e 



(17) 



In building an upwind code to solve a hyperbolic system of conservation equa- 
tions, the eigen-structure (eigenvalues and eigenvectors) of the Jacobian ma- 
trix is required. Eigen-structures for relativistic hydrodynami c s in m ultidimen- 
sions were previously described, for instance, in lDonat et al. (Il998l) . Howe ver, 
the state vector in this paper is different from that of Donat et al. ( 1998| ). so 
the eigen-structure is different. In the following, our eigen-structure of equa- 
tion (16) is presented. We first define the specific enthalpy, h, and the the 
sound speed, c s , respectively as 



e + p 



TP 
ph 



Then the eigenvalues of A x for j = x are 



(1 - Cl) V X - V /(1-^) C 2 [ 1_ U 2 C 2_ (1 _ C 2 )V 2 ] 



1 — V 2 C? 



a 2 = v x , 
Cl3 = v x , 



(19) 

(20) 
(21) 
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a 4 = v x , 



a 5 = 



(1 - c 2 ) v x + y/(l - v 2 ) c 2 [1 - v 2 c 2 s - (1 - c 2 ) t; 
1 — v 2 c 2 



(22) 
(23) 



The eigenvalues ai_5 represent the five characteristic speeds associated with 
two sound wave modes (01,5) and three entropy modes (02-4). 



The complete set of the corresponding right eigenvectors (A X R = an) is 

T 



Ri 



R9 = 



Ra 



1 ~ VxOi (1 - v x ai)v y (1 - v x ai)v z 

rh(i-vir au i-vi ' 



l-vl 



-n*h-i)v y 

, W, J-, w, u 



r 2 (2/i - 1) (t; 2 - t; 2 ) + h 



1 T 



Th 

-T(2h-l)v z 
h 



,v x ,0,0, 1 



,0,0,1,0 



1 - ^gfl5 (1 - v x a<o)v y (1 - v x a 5 )vz 

L /1 _.9\ ' °55 1 9 > 1 ..9 ' 1 



^ (1 - vi) 



l-vl 



l-vl 



(24) 
(25) 
(26) 
(27) 
(28) 



The complete set of the left eigenvectors (LA X = aL), which are orthonormal 
to the right eigenvectors, is 



-Th(v x -a 5 ) -T 2 (2h- l)(v x -a 5 )v v -T 2 (2h - 1) (v x - a 5 ) 



(h - 1) (a 1 - a 5 ) 



'12, 



(h - 1) (ai - a 5 ) 



(/i - 1) (01 - a 5 ) 



^15 



,(29) 



r/w„ [r 2 (2/i-l) (t; 2 - t; 2 ) + ft] ^ r 2 (2fe - 1) v 2 y T 2 {2h-l)v y v z 
h-V (h-l)(l-v 2 x ) ' h-1 + ' /i-l 



[r 2 (2h - 1) (t; 2 - v 2 x ) + fc] Wj 



(30) 



(h-l)(l-v 2 x ) 

Th [T 2 {2h-l){v 2 -v 2 x ) + l]v x T 2 {2h-l)v y T 2 {2h-T)v 



h-V (h-l)(l-v 

-r 2 (2h - 1) (v 2 -v 2 x )-i 

(h - 1) (1 - vl) 



h-1 



h-1 



(31) 
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Thv z [T 2 (2h - 1) (v 2 - v 2 x ) + h] v x v z Y 2 (2h - 1) v y v z T 2 (2h - 1) v\ 



-[T 2 {2h-l){v 2 -v 2 x ) + h] v z 



h-1 



(h-l)(l-v 2 ) 

—Th {y x — ai) 
(/i - 1) (a 5 - ai) 



(32) 



+ 1, 



T 2 (2/i - 1) (u, - a x ) ^ -r 2 (2h - 1) (u s - Gl ) ^ 



(/i - 1) (a 5 - ai) 



(A - 1) (a 5 - ai) 



A 



55 



(33) 



where the auxiliary variables are defined as 

A = -F 2 (2h-l)(v 2 -v 2 x ) + l} (v x -a 5 )v x + 1 
12 (h - 1) (1 - vl) {a x - a 5 ) a x - a 5 ' 

a [r 2 (2h - 1) (i; 2 - v 2 x ) + 1] (^ - gg) a 5 

[h — 1) (1 — ?;~J (ai — a 5 ) a\ — a 5 



V52 



[r 2 (2h - 1) (t; 2 - vl) + 1] („. - ai ) «, + 1 (3g) 



(h - 1) (1 - vl) (a 5 - ai) a 5 -ai' 



Ass = [r 2 (2fc - 1) (t; 2 - t; 2 ) + 1] (v x - aQ _ a_x (3?) 
(/i - 1) (1 - vl) (a 5 - ax) a 5 -ax' 



The eigenvalues and eigenvectors of A y and A 2 can be obtained by properly 
redefining indices. We note that the eigenvalues are same regardless of the 
choice of state or parameter vectors. But the right and left eigenvectors are 
different or can be presented in different forms. 



2.3 One-Dimensional Functioning Code Based on the TVD Scheme 



The TVD scheme we employ to build one-d ime nsional functionin g code is 
practically identical to that in lHarten and iRvu et al. But for 

completeness, the procedure is shown here. The state vector at the cell 
center i at the time step n is updated by calculating the modified flux vector 
fx,i±i/2 along the x-direction at the cell interface i ± 1/2 as follows: 



L x q, 
fx,i+l/2 



At n 
Ax 



/ 



x,i+l/2 



f 



x,i-l/2 J i 



fm)+fm +1 : 



Ax 5 



,At n 



Pk,i+l/2 — Qk{ A — °fc,i+l/2 + 7fe,i+l/2)"fc,j+l/2 — \9k,i + 9k,i+l) 



(38) 
(39) 
(40) 
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7M+1/2 



(9k,i+i - 9k,i) / a k,i+i/2 for a k)i+l/2 ^ 0, 
for ajfc )i+ i/ 2 = 0, 



flf fc)i = sign(^ ii+1/2 )max{0,mm[|5( M+ i / 2|,sigri(5( fcii+ i / 2)fi'/ Ci i-i/2]}, 



9k,i+l/2 



At n 

Q ' " n 



Ax k ^ 2} \ Ax k > t+1 / 



At n 



&k,i+l/2, 



Qk{x) 



+1/2 - ^k,i+l/2 ' ~ 1i ) 



x 2 /Ae k + e k for \x\ < 2e k , 
\x\ for \x\ > 2e k . 



(41) 

(42) 
(43) 
(44) 

(45) 



Here, k — 1 to 5 stand for the five characteristic modes. The internal parame- 
ters s k s are associated with numerical viscosity, and defined for < e k < 0.5; 
£^5 = 0.1 — 0.3 for the sound wave modes and e 2 -4 = — 0.1 for the entropy 
modes are reasonable choices. 

We note that the flux limiter in equation (42) is the min-mod limiter. The min- 
mod limiter is known to be very stable but has the cost of additional diffusion. 
To reproduce sharper structures with less diffusion, other flux limiters, such 
as the monotonized central difference limiter (MC limiter) 

g k>i = sign(^ ii+1/2 )max{0,min[^(|^ fcii+ i / 2| + sign(g kii+1/2 )g k ,i-i/2), 2\g k>i+1 / 2 \ 



2sign(g kti+1/2 )g kt i- 1 / 2 ]}, 



(46) 



or the superbee limiter 

g k ,i = sign(£fc ji+1 / 2 )max{0,mn[|^ 

sign(# M+1/ 2)<7 M _i/2]}, (47) 



may be used; however, these limiters are more susceptible to oscillations at 
discontinuities. In the tests described in §3, the min-mod limiter was used. 

In order to define the physical quantities at the c ell interfaces, t he TVD scheme 
originally used the Roe's linearizion technique ( Harten . 1983). Although it is 
possible to implement this linearizi on technique in the relativist ic domain in 
a computationally feasible way (see Eulderink fc Mellema . 1995| ). there is un- 
likely to be a significant advantage from the computational point of view. 
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Instead, we simply calculate the algebraic averages of quantities at two adja- 
cent cell centers to define the physical quantities at the cell interfaces; 



v x,i + v x,i+l _ Vy^i + 
V x ,i+l/2 — ^ , Vy ji+ x/2 — , f 2 ,i+l/2 



V z ,i + "^2,1+1 



(46 



l i+l/2 



C s ,i+l/2 



hi + h i+ i 



(7 - 1) ih i+1/2 - 1 



h 



i+1/2 



1/2 



(49) 
(50) 



2.4 Multidimensional Extension 



To extend the one- dimensional code to multidimensions, the procedure de- 
scribed in the previous subsection is applied separately to the y and z-directions. 
Multiple spatial dimen sions are treated through the Strang-type dimensional 
splitting ( Strang , 1968| ). Then, the state vector is updated by 



L 7 L,.L a 



(51) 



In order to maintain second-order accuracy in time, the order of the dimen- 
sional splitting is permuted as follows 

L z LyL X i L x LyL Zl L x L z Ly^ LyL z L x ^ LyL x L Zl L z L x Ly. (52) 



The time step At n is restricted by the usual Courant stability condition 

CciourA^ CcourAy CcomAz 



Af 



mm 



max Ki+i/ 2 y max Ki+i/ 2 y max «i+i/ 2 ) 



(53) 



The Courant constant should be Cc OU r < 1- We typically use Cc OU r ^ 0.9. The 
time step is calculated at the beginning of a permutation sequence and used 
through the complete sequence. 



2.5 Lorentz Transformation 



In the code, the conserved quantities D, M i} and E in the reference frame 
are evolved in time, but the physical quantities p, Vj, and e in the local rest 
frame are needed for fluxes to be estimated. The quantities p, Vj, and e can 
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be o btained through Lorentz transformation of equations (5)-(7) at each time 
step. Schneider et al. ( 1993f ) showed that the transformation is reduced to a 
single quartic equation for v 



f( v ) = [ 7U (E - Mv) — M (l — v 2 )] 2 - (l - v 2 ) v 2 (7 



l) 2 D 2 



0, (54) 



where M 2 = M 2 + M 2 + M 2 . They also showed that the physically meaningful 
solution for v is between the lower limit, i>i, and the upper limit, t>2, 



jE- ^^Ef- 4(7-1) M 
Vl = 2( 7 -l)M ' V2 = ¥' (55) 



and that the solution is unique. Once v is known, the quantities p, Vj, and e 
can be straightforwardly calculated from the following relations 

D , , 

p = y » ( 56 ) 

■U, M y M 2 

M M M y ' 

e = E — M x v x - M y v y - M z v z . (58) 



Equation (54) could be solved using a numeric al procedure such as the Newton- 
Raphson root-finding method, as suggested in lSchneider et al. (|l993l ). A prob- 
lem with the numerical approach is, however, that iterations can fail to con- 
verge. For instance, convergence can fail if one of the relativistic conditions 
is violated due to numerical errors, e.g., M > E, in a cell. This occurs 
mostly in extreme regimes. In addition, we found that convergence is often 
slow or sometimes fails in the limit M E. On the other hand, quartic 
equations have analytic s olutions. The general form of roots can be found 
in standard books such as Abramow itz fc Stegun (|l972h or on webs such as 



http://mathworld.wolfram.com/QuarticEquation.html'. Although it is too 



complicated to prove analytically, we found numerically that for the physical 
meaningful values of v and c s , v < 1 and c s < y/j — 1, among the four roots of 
equation (54), two are complex and the other two are real. While the smaller 
real root is smaller than the lower limit v±, the larger real root is between 
the two limits V\ and t>2. So the larger real root is the one we are looking for, 
and we use its analytic formula in our code. The advantages of the analytic 
approach are obvious. It always gives a solution we are looking for, and it is 
easier to predict and deal with unphysical situations if one of the relativistic 
conditions is violated due to numerical errors. 
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3 Numerical Tests 



3.1 Relativistic Shock Tube 



We have performed two sets of relativistic shock tube tests in the one, two, 
and three-dimensional computational boxes with x — [0,1], y — [0,1], and 
z = [0, 1]. Initially two different physical states are set up perpendicular to the 
direction along which waves propagate; along the x-axis in the one-dimensional 
calculation, along the diagonal line connecting (0, 0) and (1, 1) in the two- 
dimensional calculation, and along the diagonal line connecting (0, 0, 0) and 
(1, 1, 1) in the three-dimensional calculation. The initial states of the first test 
are 



(P, V X ,Vy,V Z ,p) 



(10,0,0,0, 13.3) < x, (x + y)/2, (x + y + z) /3 < 1/'. 
(1,0,0,0,10" 6 ) 1/2 <x, {x + y)/2, (x + y + z) /3 < : 



The initial states of the second test are 

(1,0,0,0,10 3 ) < x, (x + y)/2, (x + y + z) /3 < 1/2 



(p, v x ,v y ,v z ,p) 



(1,0,0,0,10- 2 ) 1/2 < x, (x + y)/2, (x + y + z) /3 < 1 



In equations (59) and (60), the expressions within inequalities are for one, two, 
and three dimensions, respectively. The first test involves a mildly relativistic 
flow and the second test involves a highly relativistic flow. In both tests, we 
assume the adiabatic index 7 = 5/3 and the outflow condition is used for 
the x, y, and ^-boundaries. Both tes ts were previously considered by several 
authors (e.g., Marti fc Miiller . Il996| ). The estimation of accuracy was done 



by compar i ng; the num erical solutions with t he exact solutions described in 
Thompson! (|l986h and iMarti fc Mimeil (|l994h . In Fi gures 1(a) and (b), our 



numerical solutions are shown as open circles and the exact solutions are 
represented by solid lines. 

Figure 1(a) shows the mildly relativistic shock tube test done using 256, 256 2 , 
and 256 3 cells with a Courant constant Cc OU r = 0.9 and the parameters £^5 = 
0.1 and £ 2 _ 4 = 0. The plots of one, two, and three-dimensions correspond to 
times t = 0.4, 0.4v^2, and 0.4a/3, respectively. Structures such as the shock 
front, contact discontinuity and rarefaction wave are accurately produced. 
There are actually slight improvements in accuracy in the multidimensional 
calculations. Figure 1(b) shows the highly relativistic shock tube test done 
again using 256, 256 2 , and 256 3 cells with a Courant constant Cc OU r = 0.6 
and the parameters £1,5 = 0.1 and £2-4 = 0. The plots of one, two, and 



three-dimensions correspond to times t = 0.4, 0.4a/2, and 0.4\/3, respectively. 
The flow is more extreme, but the structure is correctly reproduced without 
spurious oscillations. But in the rest mass density profile the peak does not 
reach the value of the exact solution due to the coarseness of computational 
cells. According to our tests, in a one- dimensional calculation, the peak can 
be accurately reproduced when 2048 numerical cells are used. There are also 
improvements in accuracy in the multidimensional calculations. 

For a more quantitative comparison, we have calculated the norm errors of the 
rest mass density, velocity, and pressure for different dimensions. The errors 
shown in Table 1 are calculated at the same times as in Figure 1. The errors 
are gradually reduced as the dimensionality increases and demonstrate a good 
agreement between the numerical and exact solutions. Note that the values of 
|| i£(p) || exceeding unity are still acceptable because these are from the initial 
large value of pressure. 



3.2 Relativistic Wall Shock 



A one-dimensional relativistic wall shock test has been performed in the com- 
putational box of x = [0, 1]. Initially a gas with extreme velocity occupying 
all numerical cells propagates along the gainst a reflecting wall placed 

at x — 1. As the gas hits the wall, it is compressed and heated and eventually 
a reverse shock is generated. The initial condition of this test is 

{p,v x ,v y ,v z ,p) = (l,0.999999, 0,0, 10~ 4 ) 0<x<l. (61) 



The adiabatic index 7 = 5/3 is assumed and the inflow boundary condition 
is use d at x = 0. It is an other test which was widely used by several authors 
(e.g.. lDonat et allll998h . 

The relativistic ju mp condition for strong sh ocks with negligible preshock 
pressure is given bv lBlandford fc McKeel (|l97fih 



P = P- 



(7 - l)Tv 

r + 1 ' 
7r + 1 



7 



0. 



(62) 

(63) 
(64) 



V 



P (r-i)( 7 r + i) 



(65) 
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Here, v s is the shock velocity and the superscript * represents the postshock 
quantities, while the quantities without any superscript refer to the preshock 
gas. 

Figure 2 shows the structure at t = 0.75 when the reverse shock is located at 
x = 0.5. The calculation has been done using 512 computational cells with a 
Courant constant Ccour = 0.9 and the parameters £1,5 = 0.3 and 52-4 = 0.1. 
The numerical solution is drawn with open circles and the exact solution is 
represented by solid lines. The numerical and exact solutions match exactly 
without any oscillation or overshoot in the rest mass density, velocity, and 
pressure profiles. 

With different inflow velocities, we have calculated the mean errors in the rest 
mass density, velocity, and pressure. The errors are calculated for the same 
time as in Figure 2 and given in Table 2. Note that the order of the mean 
errors is 10~ 3 , and that the accuracy does not depend systematically on the 
investigated Lorentz factor. The mean error in the rest mass density is < 0.5% 
for all the Lorentz factors and about 0.25% for the maximum Lorentz factor. 
This accuracy is comparable to or better than that of other published upwind 
scheme codes. 



3.3 Relativistic Blast Wave 

The propagation of a relativistic blast wave has been tested in the two- 
dimensional computational box with x = [0, 1] and y — [0, 1]. A gas of high 
density and pressure is initially confined in a spherical region and the sub- 
sequent explosion is allowed to evolve. This makes a spherical blast wave 
propagate outward. The initial condition of this test is 

|(10,0,0,0,10 3 )0<v^T^<l/2, 
[p,v x ,Vy,v z ,p) = < (66) 
[(1,0,0,0,1) outside. 

The adiabatic index is taken to be 7 = 4/3 and the reflecting and outflow 
boundary conditions are used. 

The calculation has been done using 512 2 cells with a Courant constant 
Ccour = 0.6 and the parameters e 1)5 = 0.1 and e 2 -4 = 0. To test the sym- 
metry properties of the code, the calculation has been stopped before a re- 
verse shock reaches the inner reflecting boundary. Figure 3 shows the profiles 
of the rest mass density, velocity, and pressure measured along the diagonal 
line connecting (0, 0) and (1,1) at t = 0.7. The spherical blast wave success- 
fully propagates to a larger radius, and we have found that all structures in it 
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preserve the initial symmetry. 



3.4 Relativistic Hawley-Zabusky Shock 



In order to test the applicability of the code to complex relativistic flows, we 
have performed a two-dimensional test simulation of the relati vistic version of 
the Ha wley-Zabusky shock. The test was originally suggested bv lHawlev fc Zabuskv 
(1989) ibr non-relativistic hydrodynamics. Almost the same physical values as 
in the original paper are used here. Initially a plane-parallel shock with a 
Mach number 1.2 propagates along the x-axis into two regions of different 
density. The regions are separated by oblique discontinuity whose inclination 
is 30° with respect to the x-axis. The density jumps three times across the 
discontinuity. The initial configuration is summarized as 



(p,V X ,Vy,V Z ,p) 



(1, 0.6, 0, 0, 0.48) < x < 1/16, < y < 1, 

(1, 0, 0, 0, 0.48) 1/16 < x < V3y + 1/4, < y < 1, (67) 

(3,0,0,0,0.48) outside. 



The adiabatic index 7 = 1.4 is used. The inflow and outflow conditions are used 
at the x-boundaries and the reflecting condition is used at the ^-boundaries. 

The simulation has been done in the two-dimensional computational box with 
x = [0, 8] and y = [0, 1] using a uniform numerical grid of 2048 x 256 cells. A 
Courant constant Ccour = 0.9 and the parameters Ei^ = 0.1 and £2-4 = were 
used. We have simulated this test until t = 20 in order to see the long term 
evolution. The passage of the planar shock through the discontinuity causes 
the Kelvin-Helmholtz instability to occur along the discontinuity and end up 
with formation of vortices. The vortices roll up, interact, and merge during 
the simulation; the detailed morphology and the number of vortices formed 
are somewhat sensitive to numerical resolution. Figure 4 shows the gray-scale 
images of the rest mass density at different times (t = 2, 11, and 20). Because 
all the structures are dragged to the right boundary as time goes on, only the 
left, middle, and right half of the computational box are shown at t = 2, 11, 
and 20, respectively. The vortices along the discontinuity are clearly formed 
and overall the morphology is similar to that of the non-relativistic simulation. 



3.5 Relativistic Extragalactic Jets 



Finally, in order to test the applicability of the code to realistic relativistic 
flows, we have simulated a two-dimensional relativistic extragalactic jet prop- 
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agating into homogeneous medium. The relativistic jet inflows with a velocity 
0.99 to the computational box of x = [0, 4] and y = [0, 1]. The jet has initially 
radius 1/8 (32 cells) and Mach number 8.76. The density ratio of the jet to 
the ambient medium is 0.1 and the pressure of the jet is in equilibrium with 
that of the ambient medium. The initial condition for jet inflow and ambient 
medium is summarized as 

;i, 0.99, 0, 0, 0.1) < x < 1/32, < y < 1/8, 
{p,v x ,v y ,v t ,p) = \ (68) 
'10,0,0,0,0.1) outside. 



The adiabatic index 7 = 4/3 is used. The inflow and outflow conditions are 
used at the ^-boundaries and the reflecting and outflow conditions are used 
at the ^-boundaries. 

The simulation has been done using a uniform numerical grid of 1024 x 256 
cells with a Courant constant Cc our = 0.3 and the parameters e\$ = 0.3 and 
£2-4 = 0.1. Figure 5 shows the gray-scale images of logarithm of the rest mass 
density, pressure, and Lorentz factor at t — 5 when the bow shock reaches 
the right boundary. We can clearly see the dominant structures of bow shock, 
working surface, contact discontinuity, and cocoon. It is clear that the internal 
structure of the relativistic jet is less complex compared to that of a non- 
relativistic jet due to the effects of high Lorentz factor. The overall morphology 
and dynamics of our simulatio n match roughly with those of previous works, 
e.g., Duncan fc Hug (1994), although the initial conditions and the plotted 



epoch are different. 



4 Summary and Discussion 



A multidimensional code for special relativistic hydrodynamics was described. 
It differs from previous codes in the follo wing aspects: 1) It is based on the 
total variation diminishing (TVD) scheme (Harten, 1983| ). which is an explicit 



Eulerian finite difference upwind scheme and an extension of the Roe scheme 
to second-order accuracy in space and time. 2) It employs a new set of con- 
served quantities, and so the paper describes a new eigen-structure for special 
relativistic hydrodynamics. 3) For the Lorentz transformation from the con- 
served quantities in the reference frame to the physical quantities in the local 
rest frame, an analytic formula is used. 

To demonstrate the performance of the code, several tests were presented, 
including relativistic shock tubes, a relativistic wall shock, a relativistic blast 
wave, the relativistic version of the Hawley-Zabusky shock, and a relativistic 
extragalactic jet. The relativistic shock tube tests showed that the code clearly 
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resolves mildly relativistic and highly relativistic shocks within 2 — 4 numer- 
ical cells, although it requires more cells for resolving contact discontinuities. 
The relativistic wall shock test showed that the code correctly captures very 
strong shocks with very high Lorentz factors. The relativistic blast wave test 
showed that blast waves propagate through ambient medium while preserving 
the symmetry. The test simulations of the relativistic version of the Hawley- 
Zabusky shock and a relativistic extragalactic jet proved the robustness and 
flexibility of the code, and that the code can be applied to studies of practical 
astrophysical problems. 

The strong points of the new code include the following: 1) Based on the TVD 
scheme, the code is simple and fast. The core routine of the TVD relativistic 
hydrodynamics is only about 300 lines long in the three-dimensional version. It 
runs only about 1.5 — 2 times slower than the non- relativistic counterpart (per 
time step). Yet, tests have shown that the code is accurate and reliable enough 
to be suited for astrophysical applications. In addition, the use of an analytic 
formula for Lorentz transformation makes the code robust, so it ran for all 
the tests we have performed without failing to converge. 2) The code has been 
built in a way to be completely parallel to the non-relativistic counterpart. 
So it can be easily understood and used, once one is familiar with the non- 
relativistic code. In addition, the techniques developed for the non-relativistic 
code such as parallelization can be imported transparently. 

Finally, the code is currently being applied for studies of relativistic jet in- 
teractions with inhomogeneous external media and turbulence of relativistic 
flows. The results will be reported in separate papers. 
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Table 1 



Norm errors for the relativistic shock tube tests 



RST 


"cell 


\\E(p)\\ 


\\E(v)\ 


1 


\\E(p)\\ 


(a) ID 


256 


1.1688E-01 


6.0952E- 


-02 


9.3517E-02 


2D 


256 2 


1.1264E-01 


6.0586E- 


-02 


9.6789E-02 


3D 


256 3 


9.1309E-02 


5.8222E- 


-02 


8.7047E-02 


(b) ID 


256 


1.7506E-01 


2.6591E- 


-02 


5.2191E+00 


2D 


256 2 


1.6375E-01 


1.9552E- 


-02 


4.3126E+00 


3D 


256 3 


1.3840E-01 


1.3533E- 


-02 


2.8773E+00 
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Table 2 

Mean errors for the relativistic wall shock tests 



V 


r 


E{p) 




E{v) 




E(p) 




0.9 


2.3 


4.7423E- 


-03 


3.1483E- 


-03 


5.8100E- 


-03 


0.99 


7.1 


3.1938E- 


-03 


2.3634E- 


-03 


2.5168E- 


-03 


0.999 


22.4 


3.1876E- 


-03 


2.6687E- 


-03 


2.5015E- 


-03 


0.9999 


70.7 


5.0532E- 


-03 


4.1790E- 


-03 


3.8529E- 


-03 


0.99999 


223.6 


2.8425E- 


-03 


2.4914E- 


-03 


2.1466E- 


-03 


0.999999 


707.1 


2.4855E- 


-03 


2.0237E- 


-03 


1.8747E- 


-03 
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Fig. 1. (a) ID, 2D, and 3D mildly relativistic shock tube tests. The calculations 
have been done with the initial states in equation (59) using 256, 256 2 , and 256 3 
cells. The numerical solutions (open circles) and the exact solutions (solid lines) are 
plotted at t = 0.4, OAy/2, and OAy/3. (b) ID, 2D, and 3D highly relativistic shock 
tube tests. The calculations have been done with the initial states in equation (60) 
using 256, 256 2 , and 256 3 cells. The numerical solutions (open circles) and the exact 
solutions (solid lines) are plotted at t = 0.4, 0.4y2, and OAy/3. 
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Fig. 2. One-dimensional relativistic wall shock test. The calculation has been done 
with the initial states in equation (61) using 512 cells. The numerical solution (open 
circles) and the exact solution (solid lines) are plotted at t = 0.75. 
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0.3 OA 0.6 O.fl j 
x. y 

Fig. 3. Two-dimensional relativistic blast wave test. The calculation has been done 
with the initial states in equation (66) using 512 2 cells. The numerical solution (open 
circles) are plotted at t = 0.7. 
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Fig. 4. Two-dimensional relativistic Hawley-Zabusky shock. The simulation has been 
carried out with the initial configurations in equation (67) using 2048 x 256 cells. 
Gray-scale images show the rest mass density at t = 2, 11, and 20 (top to bottom), 
using linear scales that range from 1.0 (black) to 6.75 (white). 
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Fig. 5. Two-dimensional relativistic extragalactic jet. The simulation has been car- 
ried out with the initial conditions in equation (68) using 1024 x 256 cells. Gray-scale 
images show logarithm of the rest mass density, pressure, and Lorentz factor (top 
to bottom) at t = 5, using logarithmic scales that range from —0.28 (black) to 1.98 
(white) for log(density), —1.0 (black) to 1.30 (white) for log(pressure), and (black) 
to 0.85 (white) for log(Lorentz factor). 
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